
function [ES2,EQ2,S2,EQU2] = surplus2(theta,B,D,P,xi,n,options,cmax,save_eqL,save_eqU)

N = sum(n);
NS = size(xi,1);

% PRI-PVEM joint surplus with joint PVEM candidate
ES2 = zeros(N,NS);
if save_eqL == 1
    EQ2 = zeros(N,5,NS);
    S2 = zeros(N,5,NS);
end
if save_eqU == 1
    EQU2 = zeros(N,5,NS);
end

for t = 1:N
    Xt = [blkdiag(kron(eye(2),[0,0,1,D(t,:)]),kron(eye(2),[0,1,D(t,:)]),...
        [0,0,1,D(t,:)]),log([P(t,1:2),sum(P(t,3:4)),sum(P(t,3:4)),P(t,5)]')];
    parfor ns = 1:NS
        if save_eqU == 1
            [eqL,~,eqU,~] = EQ_bounds(theta,B,Xt,xi(ns,:)',options,cmax,2); % largest and smallest spending equilibria
        else
            eqL = EQ_bounds(theta,B,Xt,xi(ns,:)',options,cmax,2);
        end
        [~,s] = Shares([B;0;0;xi(ns,:)'],zeros(5,1),[Xt,eqL,eqL.^2],0,ones(5,5),n,zeros(1,2),1);
        ESt(ns) = (theta(3)+theta(4))*log(s(3))-eqL(3);
        if save_eqL == 1
            EQt(:,ns) = eqL;
            St(:,ns) = s;
        end
        if save_eqU == 1
            EQUt(:,ns) = eqU;
        end
    end
    ES2(t,:) = ESt;
    if save_eqL == 1
        EQ2(t,:,:) = EQt;
        S2(t,:,:) = St;
    end
    if save_eqU == 1
        EQU2(t,:,:) = EQUt;
    end
end



